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In equilibrium, colloidal suspensions governed by short-range attractive and long-range repulsive 
interactions form thermodynamically stable clusters. Using Brownian dynamics computer simula¬ 
tions, we investigate how this equilibrium clustering is affected when such particles are self-propelled. 
We hnd that the clustering process is stable under self-propulsion. For the range of interaction 
parameters studied and at low particle density, the cluster size increases with the speed of self- 
propulsion (activity) and for higher activity the cluster size decreases, showing a non-monotonic 
variation of cluster size with activity. This clustering behaviour is distinct from the pure kinetic 
(or motility-induced) clustering of self-propelling particles which is observed at significantly higher 
activities and densities. We present an equilibrium model incorporating the effect of activity as 
activity-induced attraction and repulsion by imposing that the strength of these interactions de¬ 
pend on activity superlinearly. The model explains the cluster size dependence of activity obtained 
from simulations semi-quantitatively. Our predictions are verifiable in experiments on interacting 
synthetic colloidal microswimmers. 

PACS numbers: 82.70.Dd, 64.75.Xc 


I. INTRODUCTION 

Concentrated colloidal or protein solutions which are 
governed by a combination of short-range attractive and 
long-range repulsive interaction potentials exhibit a sta¬ 
ble clustering phenomenon in equilibrium at finite tem¬ 
perature and moderate densities. This phenomenon is 
first predicted by theory [T] and simulation ma and has 
been confirmed in experiments laiT]. The intuitive expla¬ 
nation [T] for equilibrium clustering lies in the fact that 
the short-ranged attraction first leads to growth of clus¬ 
ters in an initially dilute suspensions of particles. The 
growth stops, however, when the cluster size reaches a 
characteristic size where the long-ranged repulsion leads 
to an increase in the self-energy of the cluster. In equi¬ 
librium, at finite temperature, this leads to a typical av¬ 
erage cluster size which depends on the interaction pa¬ 
rameters and the imposed global particle density. While 
the details of this equilibrium cluster process are under¬ 
stood for a decade by now, recent developments have 
considered self-propelled (or active) particles which dis¬ 
sipate energy leading to synthetic microswimmers 019]. 
These particles also exhibit a purely kinetic clustering if 
the strength of self-propulsion is sufficiently large 01101 
which has recently been found in experiments [IIHIOI and 
explored by simulation [T2j [T4]-[22] and theory pn[2Q| l23l- 
[28]. This purely motility-induced clustering occurs for 
repulsive systems and is therefore absent in equilibrium 
(i.e. for vanishing drive). The study of Redner et al [29] 
showed reentrant phase behavior in active Lennard-Jones 
particles, where attractive interaction and activity com- 
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pete to stabilize phase-separated states at low and high 
activities, respectively. 

Here we link the fields of equilibrium clustering to that 
of microswimmers. We consider the equilibrium cluster¬ 
ing and study how this is affected by an imposed self¬ 
propulsion. The motivation to do so is three-fold. First 
of all, this is an interesting problem in itself since upon in¬ 
creasing the self-propulsion, there are two counterbalanc¬ 
ing effects: on the one hand, the self-propulsion leads to a 
higher mobility and hence an effect which is expected to 
correspond to an increase of temperature. On the other 
hand, however, the self-propulsion yields a larger sticking 
probability of neighbouring particles which would favor 
and enhance the clustering tendency. The second moti¬ 
vation comes from the fact that one needs to understand 
whether there is a hidden pathway between the two dif¬ 
ferent kinds of clustering mentioned above, i.e. to check 
whether they are distinct or interconnected in a certain 
parameter space. Finally, artificial colloidal model mi¬ 
croswimmers can be prepared with controlled interac¬ 
tions, e.g. by adding depletants m, tuned van der Waals 
attractions or charging the particles such that model col¬ 
loidal swimmers can be prepared, in principle, with short- 
ranged attraction and long-ranged repulsion. The addi¬ 
tional tunability of the interparticle potential then allows 
to control the degree of clustering which needs a system¬ 
atic understanding. 

In this paper we simulate a two-dimensional model 
of microswimmers with competing interactions by us¬ 
ing Brownian dynamics computer simulations. We use a 
model proposed by Sear et al m and Imperio and Reatto 
[3| , for which the equilibrium clustering behaviour is well- 
understood in two-dimensions but supplement this here 
for an additional self-propulsion in the simplest form by 
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neglecting explicit alignment and hydrodynamic interac¬ 
tions [32]. As a result, we find indeed that the trends 
of clustering depend on the interaction parameters. The 
self-propulsion can either increase or decrease the cluster 
size. In fact there is a complex and maybe unanticipated 
non-monotonic behaviour of the cluster size as a func¬ 
tion of increasing self-propulsion: it can first increase and 
then decrease again. This cannot be understood by sim¬ 
ple temperature rescaling as has been previously noted 
for active systems in the context of freezing [32] as well 
as in a trapping [33] and in a gravitational field m- 
Dynamic (or purely motility-induced) clustering also oc¬ 
curs in our model although at much larger drives, where 
the details of the interactions become irrelevant. In this 
case, the cluster sizes are rather small compared to that 
of equilibrium clusters. Thus the two clustering phenom¬ 
ena appear quite distinct. 

II. MODEL AND SIMULATION 



ria 

FIG. 1. (Colour online) Overall interaction potential u(r) 
showing competing attractive and repulsive interactions for 
Ca = 25k bT from Eq. ^ The inset shows the potential near 
contact. 


We model short-range attractive and long-range repul¬ 
sive interactions using a modified Lennard-Jones poten¬ 
tial ui{r) and a double-exponential potential U 2 {r) that 
was introduced previously to explain the formation of 
finite sized clusters and stripes of nanoparticles at air- 
water interface [31]. The overall potential is given as: 

u{r) = ui{r) + U 2 {r) (1) 


Here ui{r) is defined as 


ui{r) = Aclj 




( 2 ) 


and 1^2 (r) is given by 
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Here r is the interparticle distance, and a is the diameter 
of the particle. €lj and aLj are the parameters of the 
modified Lennard-Jones potential. We fix the potential 
parameters to = cr, Ra = a^ Rr = 2a and = 
e^. In the following, we use dimensionless quantities and 
express energy in units of ksT^ length in units of a, time 
in units of r = a^ jD. Here, is Boltzmann constant, 
T is temperature and D is the diffusion coefficient of a 
single passive particle. We further fix e^j = 0.0025k bT. 
Figshows the variation of potential with interparticle 
distance for = 25k bT. Note that the repulsive part of 
the double-exponential potential is rather long ranged. 

Brownian dynamics simulations in two dimensions in 
the x^-plane are performed with particles interacting via 
the potential given in Eq. 0. We simulate N = 1050 
particles using a square box with periodic boundary con¬ 
ditions. To mimic self-propulsion, the particles are de¬ 
fined with an orientation diffusing freely about the 
perpendicular 2 : axis with rotational diffusivity In 


two dimensions, the components of are given as = 
(cossin(/9^). In addition to translational Brownian 
motion, the particles are driven with constant speed v 
along their orientation e^. Here there are no aligning 
interactions as the pair potential is independent of ori¬ 
entations. Moreover hydrodynamics interactions are ne¬ 
glected. The resulting equations of motion for the par¬ 
ticle positions {r^} and orientations {ef} are then given 
by 
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‘Pi = ( 5 ) 

Here, U = “ r^j) is the total pair potential. 

The self-propulsion speed of the particle is referred to in 
terms of the dimensionless Peclet number Pe defined as 


The Gaussian noise models the stochastic solvent 
kicks. It has a zero mean and variance {t')) = 

25ijDl6{t — t'), where 1 is the identity matrix. Simi¬ 
larly, the stochastic random torque has a zero mean 
and a variance of {^1 {t)^j {t')) = 2DrSij6{t — t'). The 
rotational diffusivity is taken as Dr = ZDja^ ^ which is 
a valid approximation for a spherical particle undergoing 
free rotational diffusion. The equations of motion Eq. 
0 and 0 are numerically integrated with a time step 
of 10“^r. The long range potential Eq. 0 is truncated 
at r = 15(7. Simulations were done for various reduced 
areal densities (typically of 0.1Za~^ unless stated other¬ 
wise) and for different values of The simulations are 
performed starting from a random initial configuration 
of particles. Typically the system is simulated for 500r 
(i.e. 5 X 10^ steps) to attain a steady-state, followed by 



















3 


production runs of another 500r. The simulations are 
replicated 5 times with different initial random configu¬ 
rations and the properties are time-averaged over all the 
replicas. Careful tests were performed to check that the 
system achieved a steady state by monitoring the satu¬ 
ration of the cluster size as a function of time and by 
making sure that there is enough exchange dynamics be¬ 
tween the clusters. 

To define a cluster we use a cutoff distance of 1.5cr, 
which corresponds to the interparticle distance where the 
potential is roughly half of the potential at contact. Two 
particles belong to the same cluster if they are connected 
by a sequence of other particles which are all separated 
by less than 1.5cr. The average cluster size is calculated 
from 


N 

<n>=^nP(n) (7) 

n=l 

Here, P{n) is the probability to find a cluster of n num¬ 
ber of particles at steady state. We also monitor the 
fluctuations in the cluster size by calculating the reduced 
variance 


Var{n) 


<in? > — <n>‘^ 
< n 


(8) 


III. CLUSTERING OF ACTIVE PARTICLES 

Passive particles with short-ranged attractive and 
long-ranged repulsive interactions defined in Eq. Q 
show equilibrium clustering at certain densities. In par¬ 
ticular, equilibrium clustering occurs for a density of 
about 0.13 and for attraction energies in the range 
of 12kBT to 25 /cbT [3]. These clusters are quasi-circular 
in shape. The effect of activity on clustering is shown 
in Fig. where the average cluster size < n > is shown 
versus the Peclet number Pe. For the cases of vanish¬ 
ing activity, the average cluster size is about < n 14 
independent of e^. 

Strikingly, the dependence of the average cluster size 
on activity is non-mo notonic. The cluster size first in¬ 
creases with increasing activity and attains a maximum 
before it finally decreases at higher activity. This trend 
is seen for all the values studied here. The critical ac¬ 
tivity corresponding to the maximum of cluster size in¬ 
creases with increasing e^. This behaviour points to the 
possibility that the activity manifests itself as an effective 
attraction which increases the cluster size until the ac¬ 
tivity gets so high that particles are eventually removed 
from the cluster overcoming attractive interactions be¬ 
tween the particles in the cluster. It is interesting that 
we find a maximum in the cluster size at intermediate 
Pe while Redner et al [29] find a suppression of phase 
separation at intermediate Pe. Therefore, our findings 
are qualitatively different from that of Redner et al [29| 
due to the combination of attraction and repulsion. 



Pe 


FIG. 2. (Colour online) Effect of activity on the average 
cluster size < n > for various values of Ca- 


Representative snapshots showing clusters for Pe = 0 
and Pe = 6 (at fixed = 2bkBT) are shown in Fig. 
and i3, respectively. The clusters exhibit an inner crys¬ 
talline structure as found in previous simulations for dy¬ 
namical clustering m- Moreover, there is a large spread 
in cluster sizes. This is also documented by the normal¬ 
ized size distribution function P(n) which is shown in the 
insets of Fig. and lb- An increased activity leads to 
a much larger distribution in the cluster size at steady- 
state. This is documented by the reduced variance of the 
cluster size distribution which steeply increases with Pe 
as shown in Fig. until it reaches a maximum and de¬ 
creases as it is correlated with the average cluster size. 


IV. EFFECTIVE EQUILIBRIUM MODEL 


We now present a phenomenological model to explain 
the clustering behavior observed in active particles by 
balancing interparticle interactions with activity. First 
we consider passive particles with short-ranged attractive 
and long-ranged repulsive interactions and include the ef¬ 
fect of activity in the model later. Consider monomeric 
discs of diameter a assembled as circular clusters of uni¬ 
form diameter d with n number of monomers. The discs 
interact with themselves via a short-range attraction and 
long-range repulsion. Following the approach of Groe- 
newold and Kegel [T] , the free energy of such a cluster of 
passive particles can be written as 


fin) 
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FIG. 3. (Colour online) Representative snapshots of clusters for Pe = 0 (a) and Pe = 6 (b) for Ca = 25 ksT. The insets show 
the cluster size distribution P{n) in the steady state. 



counts for attractive energy assuming the attraction to 
be short ranged and therefore scales linear in n. The 
last term is energy due to line tension A of the cluster 
boundary. Also note that n and d are related via, 



( 10 ) 


where a is the cross sectional area of the monomer (a = 
7rcr^/4). Combining Eqs. and ( [Tq| ) we get 


f{n) 

n 


pji _ 1 

X -gEr^/n - EaP 2Xx/Ea—= (11) 
V a 


Minimizing Eq. @ with respect to n yields an equilib¬ 
rium cluster size n* 


n 
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FIG. 4. (Colour online) Effect of activity on the reduced 
variance of cluster size for various ta values. 


Approximating the line tension as A « Eajo^ the equi¬ 
librium cluster size is given as 


2gEr 


(13) 


where ^ is a geometric parameter related to circular 
shape, Er and Ea are typical repulsive and attractive 
energies of a particle inside the cluster. The first term 
accounts for total repulsive energy and is of order due 
to long-range nature of repulsion, the second term ac- 


Eq. (13) gives the effect of interaction parameters 
on the equilibrium cluster size of passive particles with 
short-ranged attractive and long-ranged repulsive inter¬ 
actions. We use Eq. (13) as a starting point to ana¬ 
lyze active particles, wherein each particle propels with 
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a speed of v and rotates freely in two dimensions in ad¬ 
dition to their short-ranged attractive and long-ranged 
repulsive interparticle interactions. The critical issue is 
to know how to incorporate the effect of activity in terms 
of effective attractive and repulsive interactions in Eq. 
(13). Consistent with previous work HEniEi] , we pro¬ 
pose that the role of activity in affecting inter-particle 
interactions has the following features: 

1. Both effective attraction and repulsion increase 
with Pe. 


2. For small Pe, the increase of an activity-induced 
effective attraction is more pronounced than the 
activity-induced effective repulsion. 

3. For large Pe, the effective activity-induced repul¬ 
sion is getting more pronounced than the activity- 
induced effective attraction. 


Therefore, we replace Ea in Eq. (13) to account for an 
effective activity-induced attraction as Ea + aPe^ with 
two fit parameters: an amplitude a and an exponent p. 
Similarly we add to Er an activity-induced repulsion, i.e. 
we replace Er by Er + hPe^ with two fit parameters, 
namely an amplitude b and an exponent q which is larger 
than p. Therefore 


71(7 {E a + aPe^) 
2g{Er + bPe^) 


(14) 


We find reasonable fits for our simulation data when fix¬ 
ing p = 2 and g = 4 adjusting only the amplitudes a and 
b. The actual fit parameters are given in Table I and in 
Fig. [^the comparison between the model and simulation 
data is shown. Good fits are obtained for small but 
deviations are visible for larger Here, the values of 
Ea and Er correspond to the minimum and maximum 
of the potential described in Eq. 0 . The parameter 
^ = 1.83 is chosen such that the model predicts the size 
of the cluster size in the absence of activity. Hence this 
phenomenological model can give account for the trends 
that at low Pe, < n > increases and then decreases for 
high Pe . 


TABLE I. Fit parameters used in the effective equilibrium 
model 


taifsT) 

EaiksT) 

EriksT) 

a 

b 

12 

2.6 

0.18 

1.0605 

0.0142 

15 

3.3 

0.21 

0.5928 

0.0030 

20 

4.4 

0.31 

0.4550 

0.0008 

25 

5.5 

0.38 

0.3226 

0.0002 


V. EQUILIBRIUM VERSUS DYNAMICAL 
CLUSTERING 

Next we discuss the state behaviour of active parti¬ 
cles with competing interactions to understand the link 



Pe 

FIG. 5. (Colour on line) Comparison between predictions of 
the model (Eq. with simulation data. 


between equilibrium clustering and kinetic clustering in¬ 
duced by activity. Representative configurations ob¬ 
tained from simulations for different density and activ¬ 
ity are shown in Fig. and a corresponding state di¬ 
agram is presented in Fig. In the absence of activ¬ 
ity (Pe=0), the particles form quasi-circular clusters at 
low density and extended worm-like clusters at higher 
densities. Upon increasing activity the phase behaviour 
changes significantly. For instance, up to a reduced den¬ 
sity of 0.2, the activity increases the cluster sizes until a 
critical activity is reached and beyond this activity, clus¬ 
ters dissolve leading to a disordered fluid phase. In the 
density range of 0.3 to 0.4, we see a similar behaviour but 
now the elongated worm-like clusters are getting smaller 
in length upon increasing the activity. Further increase 
in activity again leads to a disordered fluid in this density 
regime. At higher densities such as 0.5, the disordered 
fluid is followed by phase-separation due to kinetic clus¬ 
tering induced by motility as found in earlier reports on 
purely repulsive particle [num. At these high Peclet 
numbers, details of the interparticle interaction except 
for the repulsive core are not relevant. The dynamic 
clustering observed at high Pe is therefore explained as a 
pure kinetic effect. It is well separated from the equilib¬ 
rium clustering considered earlier which occurs at small 
Pe and small densities demonstrating that these two clus¬ 
tering effects are qualitatively distinct. 

We now revisit the cluster size distribution function 
P{n) in the limit of high Peclet numbers where only the 
repulsive core is relevant, which is then to be compared to 
the previous data shown in Fig. Results for Pe = 100 
are shown in Fig. on a semi-log scale. The average 
cluster size in this case is 1.8. The clusters are distinctly 
different, both in size and shape, from that observed in 
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FIG. 6. (Colour online) Representative configurations of various states of active particles with competing interactions for varied 
Peclet number Pe and reduced number density pa^ for Ca = 25k bT 
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FIG. 7. (Colour online) State diagram of active particles with 
competing interactions in the plane spanned by Peclet number 
Pe and reduced density pa^ for Ca = 25kBT. 


equilibrium (Fig. [^) and at moderate activity (Fig. [^), 
although Ca and density are the same. This means that 
activity can be used as a knob, to an extent, in tuning 
the cluster size either to increase or to decrease depending 
upon the fixed parameters of interactions. This unique 
feature is not present for purely motility-induced clus¬ 
tering in repulsive colloids. In more detail, as can be 
deduced from Fig. the cluster size distribution P(n) 
is almost linear in the semi-logarithmic plot, except for 
small n (n < 4). This shows that there is an exponential 
decay in P{n) for large n. The data for n < 5 are com¬ 
patible with a power-law scaling. Majumdar et al. [35] 


showed that for an effective aggregating and fragment¬ 
ing particle system P{n) decays exponentially for low 
aggregation rates and decays as power-law for high ag¬ 
gregation rates. In the present case, the initial power-law 
decay denoting an effective aggregation process indicates 
motility-induced kinetic clustering at high Peclet num¬ 
bers. The exponential decay denoting effective fragmen¬ 
tation indicates fragmentation due to activity. Therefore, 
high activity induces both aggregating and fragmenting 
processes. In contrast, in the case of Pe = 6, where 
both the interaction potential and activity affect cluster¬ 
ing, the cluster size distribution P{n) does not show any 
conclusive scaling. 

Furthermore we note that our results are different from 
flying crystals which were found to exist with and without 
cohesive forces [361438] . In these studies, aligning forces 
between the particles are relevant. In our simulations, we 
have not considered aligning forces and therefore we have 
not observed flying crystals. The velocity vectors of the 
particles in the clusters obtained from our simulations are 
randomly oriented as the particles are freely rotating with 
their rotational diffusion coefficient. Therefore, there is 
only a random and undirected migration of the clusters. 


VI. CONCLUSION 

In conclusion, by using Brownian dynamics computer 
simulations, we explored how equilibrium clustering is 
affected for self-propelled colloidal particles. While this 
clustering process is stable under self-propulsion, depend¬ 
ing on the values of interaction parameter, the cluster size 
can initially increase with the strength of self-propulsion, 
before it decreases for large activity. This allows to con- 
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cle interactions. A phenomenological model is shown to 
qualitatively explain the non-monotonic variation of clus¬ 
ter size with activity. For the future, it would be inter¬ 
esting to construct a dynamical density functional theory 
for the clustering considered here by unifying the den¬ 
sity functional theory designed for equilibrium clustering 
j39] with that designed for kinetic clustering HZ]. More¬ 
over, hydrodynamic effects should be explored by more 
sophisticated simulation models m- Furthermore our 
predictions can in principle be verified in experiments on 
synthetic or bacterial microswimmers with well-defined 
interactions. In particular the combination of depletants, 
particle charge and magnetic dipole moments [m gg 
opens new ways to steer the interparticle interactions be¬ 
tween swimmers and therefore the details of the cluster¬ 
ing behaviour. 
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